#ifndef _MBPT_H_
#define _MBPT_H_
#include"help.h"
#include"twobodystate.h"
#include"system_hf.h"
#include"system_gammow.h"
#include<cmath>
#include<eigen3/Eigen/Eigen>//matrix library


using std::pow;

using Eigen::Matrix;
using Eigen::Dynamic;
using Eigen::MatrixXcd;
using Eigen::SelfAdjointEigenSolver;
using Eigen::ComplexEigenSolver;


template<typename Scalar>
void quickSort(Matrix<Scalar,Dynamic,1> & eigvals,int startPos, int endPos) 
{ 
  int i,j; 
  Scalar key; 
  key=eigvals[startPos];
  i=startPos; 
  j=endPos;
  while(i<j) 
    { 
      while(!(eigvals[j]<key)&& i<j) --j; 
      eigvals[i]=eigvals[j];
      while(eigvals[i]<key && i<j) ++i; 
      eigvals[j]=eigvals[i];
    } 
  eigvals[i]=key;
  if(i-1>startPos) quickSort(eigvals,startPos,i-1); 
  if(endPos>i+1) quickSort(eigvals,i+1,endPos); 
} 

  //sort algorithm, sort according to eigvals and reorder the eigvecs respectively
template<typename Scalar>
void quickSort(Matrix<Scalar,Dynamic,1> & eigvals,Matrix<Scalar,Dynamic,Dynamic> &eigvecs,int startPos, int endPos) 
{ 
  int i,j; 
  Scalar key; 
  key=eigvals[startPos];
  Matrix<Scalar,Dynamic,1> Vec=eigvecs.col(startPos);
  i=startPos; 
  j=endPos;
  while(i<j) 
    { 
      while(!(eigvals[j]<key)&& i<j) --j; 
      eigvals[i]=eigvals[j];
      eigvecs.col(i)=eigvecs.col(j);
      while(eigvals[i]<key && i<j) ++i; 
      eigvals[j]=eigvals[i];
      eigvecs.col(j)=eigvecs.col(i);
    } 
  eigvals[i]=key;
  eigvecs.col(i)=Vec; 
  if(i-1>startPos) quickSort(eigvals,eigvecs,startPos,i-1); 
  if(endPos>i+1) quickSort(eigvals,eigvecs,i+1,endPos); 
}



struct Channel
{
  typedef complex<double> DataType;
  //  typedef double DataType;
  Channel(){}
  Channel(int _Tz,int _Par,int _J):Tz(_Tz),Par(_Par),J(_J){}
  int Tz,Par,J;
  vector<int> bra;
  vector<int> ket;
  Matrix<DataType,Dynamic,Dynamic> Mat;
};


//perturbation for closed shell and effective interaction for open shell//
class MBPT
{
 public:
  typedef complex<double> DataType;
  //  typedef double DataType;
  typedef System_Tz_HO SystemType;
  //typedef System_Tz_HF SystemHFType;
  typedef System_Tz_Gammow SystemHFType;  
  typedef TwoBodyState_TzParJ TwoBodyStateType;
  typedef SystemHFType::OrbitalType HFOrbital;
  typedef vector<HFOrbital> HFOrbitals;

  typedef Matrix<DataType,Dynamic,Dynamic> Mat;
  typedef Matrix<DataType,Dynamic,1> Vec;
  //this constructor is only for closed shell nuclei MBPT cal. , take interaction under HF. basis.
 MBPT(const SystemHFType& _HFSystem):HFSystem(_HFSystem),HFOrbs(HFSystem.Orbitals),beta(0),realA(_HFSystem.realA)
    {
      totalOrbitals=HFOrbs.size();
      FermiSurface=HFSystem.FermiSurface;
      ModelSpace.clear();
      ModelSpace.resize(totalOrbitals,false);
      IndexInP.clear();
      IndexInP.resize(totalOrbitals,-1);      
      setup2BStates();
      // setupVMat();

      beta=0;
    }

  //read norm-ordered Hamiltonion
 MBPT(const char * orbfile,const char *intfile):HFOrbs(HFSystem.Orbitals)
    {
      HFSystem.setupFromFile(orbfile,intfile);
      realA=HFSystem.realA;
      totalOrbitals=HFOrbs.size();
      FermiSurface=HFSystem.FermiSurface;
      ModelSpace.clear();
      ModelSpace.resize(totalOrbitals,false);
      IndexInP.clear();
      IndexInP.resize(totalOrbitals,-1);
      setup2BStates();
      setupVMat();
      beta=0;
    }

  
 MBPT(SystemType*Ptr,int _A,int _Z):pSystem(Ptr),HFSystem(Ptr,_A,_Z),HFOrbs(HFSystem.Orbitals),totalOrbitals(Ptr->Orbitals.size()),ModelSpace(Ptr->Orbitals.size(),false),IndexInP(Ptr->Orbitals.size(),-1),A(_A),Z(_Z),beta(0),realA(_A){}
  void setup(int _realA=-1)
  {
    if(_realA==-1)
      realA=A;
    else
      realA=_realA;
    HFSystem.setAZ(A,Z,realA);
    HFSystem.setbeta(beta);
    HFSystem.setup();
    totalOrbitals=HFOrbs.size();
    FermiSurface=HFSystem.FermiSurface;
    CoreSumSPE=0;
    for(int i=0;i<=FermiSurface;i++)
      CoreSumSPE+=HFOrbs[i].e;
    setupPGroups();
    setupP2BStates();
    
    setup2BStates();
    //    setupVMat();
  }
  DataType get2B(int a,int b,int c,int d,int J) const
  {
    SystemHFType::TwoBodyMatEle MatEle;    
    //HFSystem.get2BmatAt(a,b,c,d,J,MatEle);
    HFSystem.get2BmatOnNeed(a,b,c,d,J,MatEle);
    if(MatEle.empty()) return 0;
    if(a==b) MatEle*=sqrt(2.0);
    if(c==d) MatEle*=sqrt(2.0);
    return MatEle[0]+(beta-1)*HFSystem.hbar_omega/realA*MatEle[3]-beta*HFSystem.hbar_omega/realA*MatEle[2];
    
    //return MatEle[0]-MatEle[3]*HFSystem.hbar_omega/realA+beta*MatEle[1]*HFSystem.hbar_omega/realA;    
  }
    
  DataType cross_coupled1(int a,int b,int c,int d,int J) const;///< type 1 cross coupled matrix
  DataType cross_coupled2(int a,int b,int c,int d,int J) const;///< type 2 corss coupled matrix
  
  void calOccuNum();
  vector<DataType> OccuNum;


  void printHFOrbitals(int num=-1) const;
  void printPOrbitals() const;
  void setPSpace(int orb)
  {
    ModelSpace[orb]=true;
    POrbitals.push_back(orb);
    IndexInP[orb]=POrbitals.size()-1;
  }
  void clearModelSpace();
  vector<bool> ModelSpace;///< size=totalOrbitals, ModelSpace[i]=true means the i th orbital belong to model space P
  vector<int> IndexInP;///< size=totalOrbitals, IndexInP[i] return the index in POrbitals of ith orbital
  vector<int> POrbitals;///< P space orbitals, POrbitals[i]  is the index in HFOrbs of ith P space orbital  
  void setupPGroups();  
  vector< vector<int> > PGroups;
  vector<int> GroupIndice;
  vector<int> OrbIndice;

  DataType CoreSumSPE;

  //solve core
  DataType diagram0B_2p2h_D1(DataType w);
  void get0BQDQ(DataType w,DataType & Q,DataType & DQ);
  DataType diagram0B(DataType w);
  DataType get0BZ(DataType w);
  DataType get0BQ(DataType w);
  void solve0B();
  DataType E0;
  DataType rhoQ0B;


  //solve onebody
  DataType diagram1B(int a,int c,DataType w);
  DataType get1BQ(int a,int c,DataType w);
  void get1BQ(int nb,Mat & Q,DataType w);
  void get1BQDQ(int nb,Mat &Q,Mat & DQ,DataType w);  
  void get1BZ(int nb,Mat & S,DataType w);
  void solve1B();
  vector<Vec> OneBodyVals;
  vector<Mat> OneBodyVecs;
  vector<Mat> P1BMat;
  vector<Vec> rhoQ1B;
  void printPGroups() const;
  void printP1BMat() const;  


  //solve twobody
  DataType diagram1B_to_2B(int a,int b,int c,int d,int J,DataType w);
  DataType diagram2B(int a,int b,int c,int d,int J,DataType w);
  DataType get2BQ(int bra,int ket,DataType w);
  void get2BQ(int nb,Mat & Q,DataType w);
  void get2BQDQ(int nb,Mat & Q,Mat & DQ,DataType w);
  void get2BZ(int nb,Mat & Z,DataType w);
  void solve2B();
  DataType mat1B_to_2B(int a,int b,int c,int d,int J);
  void mat2Bminus1B();
  void symmetrize();//similarity transform to symmetrize the final effective interaction.
  vector<Vec> TwoBodyVals;
  vector<Mat> TwoBodyVecs;
  vector<Vec> rhoQ2B;
  void printP2BMat() const;

  void printsp(const char *spfile);
  void printint(const char *intfile);
  void printQ(const char*qfile);
  
  void setupP2BStates();
  void printP2BStates() const;
  void printP2BChannels() const;
  vector<TwoBodyStateType> P2BStates;// P space 2B states;
  vector< vector<int> > P2BChannels;
  vector<Mat> P2BMat;

  void printResult() const;
  //graph plot output
  void print0BGraphData(ostream & fout,double start,double step,int n);
  void print1BGraphData(ostream & fout,int nb,double start,double step,int n);
  void print2BGraphData(ostream & fout,int nb,double start,double step,int n);

  void setbeta(double _beta){beta=_beta;}
  void debug()
  {
  }
  
  const SystemType * pSystem;
  SystemHFType HFSystem;
  const HFOrbitals & HFOrbs;
  int FermiSurface;
  int totalOrbitals;

  double beta;
  int A,Z;//core A and Z;
  int realA;
  // private:
  // below is to save matrix element before diagram cal.

  
  DataType getphmhpm(int h1,int p1,int h2,int p2,int J) const;
  DataType getphmphm(int h1,int p1,int h2,int p2,int J) const;


  void setup2BStates();
  vector<TwoBodyStateType> TwoBodyStates_hh;
  vector<TwoBodyStateType> TwoBodyStates_pp;
  vector<TwoBodyStateType> TwoBodyStates_hpm;  
  vector<TwoBodyStateType> TwoBodyStates_phm;

  void setIndice_hh(const TwoBodyStateType & State_hh,int channel,int indice_hh);
  void setIndice_pp(const TwoBodyStateType & State_pp,int channel,int indice_pp);

  vector< vector<int> > ChannelIndices_hhpp;
  vector< vector<int> > StateIndices_hhpp;//(mapab_e(a,b),J-Jmin) return the channel and pos index in V_hhpp of twobody state (ab,J)

  void setIndice_hpm(const TwoBodyStateType & State_hpm,int channel,int indice_hpm);
  void setIndice_phm(const TwoBodyStateType & State_phm,int channel,int indice_phm);

  vector< vector< vector<int> > > ChannelIndices_hpm;
  vector< vector< vector<int> > > StateIndices_hpm;
  
  vector< vector< vector<int> > > ChannelIndices_phm;
  vector< vector< vector<int> > > StateIndices_phm;
  
  void setupVMat();
  vector<Channel> V_hhhh, V_hhpp,V_pppp,e_hhpp,V_phmhpm,V_phmphm,e_phmhpm;

  int block2B;

  //a lot of diagrams

  //------------zero body 
  DataType diagram0B_2p2h(DataType w);
  DataType diagram0B_4p2h(DataType w);
  DataType diagram0B_2p4h(DataType w);
  DataType diagram0B_3p3h(DataType w);
  //----------

  
  //------------one body second order
  DataType diagram1B_2p1h(int a,int c,DataType w);
  DataType diagram1B_3p2h(int a,int c,DataType w);
  //------------

  
  //------------twobody second order
  DataType core_polarization(int a,int b,int c,int d,int J,DataType w);
  DataType diagram2B_3p1h(int a,int b,int c,int d,int J,DataType w);
  DataType diagram2B_2p(int a,int b,int c,int d,int J,DataType w);
  DataType diagram2B_4p2h(int a,int b,int c,int d,int J,DataType w);  
  //------------
  
};
#endif
